COMMENT ⊗   VALID 00002 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	 This extracts some useful information about the floating point
C00012 ENDMK
C⊗;
;;; This extracts some useful information about the floating point
;;; hardware for the tests later:
;;;	*ibeta*		radix of the floating-point representation
;;;	*it*		number of base *ibeta* bits in the significand
;;;	*irnd*		0 if floating-point addition chops
;;;	         	1 if floating-point addition rounds
;;; 	*ngrd*		number of guard bits digits for multiplication
;;;	*machep*	the largest negative integer such that
;;;			1.0+float(*ibeta*)↑*machep* is not = to 1.0
;;;	*negep*		The largest negative integer such that
;;;			1.0-float(*ibeta*↑*negep* is not = to 1.0
;;;	*iexp*		the number of bits (decimall places if *ibeta* = 10)
;;;			reserved for the exponent
;;;	*minexp*	the largest in magnitude negative integer such that
;;;			float(*ibeta*)↑*minexp* is a positive floating-point
;;;			number
;;;	*maxexp*	the largest positive integer exponent for
;;;			a finite floating-point number
;;;	*eps*		the smallest positive floating-point number such that
;;;			1.0+*eps* is not = to 1.0
;;;	*epsneg*	a small positive floating-point number such that
;;;			1.0-*epsneg* is not = to 1.0
;;;	*xmin*		the smallest non-vanishing floating-point power
;;;			of the radix
;;;	*xmax*		the largest finite floating-point number.
;;; See Cody and Waite, Software Manual for the Elementary Functions, for
;;; details.

(declare (special *ibeta* *it* *irnd* *ngrd* *machep* *epsneg* *maxexp*
		  *negep* *iexp* *minexp* *maxeps* *eps* *xmin* *xmax*))

(declare (flonum *xmax* *xmin* *epsneg*)
	 (fixnum *ibeta* *it* *irnd* *machep* *minexp* *iexp* 
		 *negep* *ngrd* *maxexp*))

(defun machar ()
 (declare (flonum a b beta betain betam1 one y z q aa bb yy xmin)
	  (fixnum i iz j k maxexp mx it kk negep))
 (let ((a 0.0)(b 0.0)(beta 0.0)(betain 0.0)(betam1 0.0)
	      (one 1.0)(y 0.0)(i 0)
	      (k 0)(mx 0))
      (do ((aa (+$ one one) (+$ aa aa)))
	  ((not (zerop (-$ (-$ (+$ aa one) aa) one)))
	   (setq a aa)))
      (do ((bb (+$ one one)(+$ bb bb)))
	  ((not (zerop (-$ (+$ a bb) a)))
	   (setq b bb)))
      (setq *ibeta* (fix (-$ (+$ a b) a)))
      (setq beta (float *ibeta*))
      (do ((it 1 (+ it 1))
	   (b (*$ one beta) (*$ b beta)))
	  ((not (zerop (-$ (-$ (+$ b one) b) one)))
	   (setq *it* it)))
      (setq *irnd* 0)
      (setq betam1 (-$ beta one))
      (cond ((not (zerop (-$ (+$ a betam1) a)))
	     (setq *irnd* 1)))
      (setq *negep* (+ *it* 3))
      (setq betain (//$ one beta))
      (do ((i *negep* (1- i))
	   (aa (*$ one betain) (*$ aa betain)))
	  ((zerop i) (setq a aa)))
      (let ((a a))
	   (do ((aa a (*$ aa beta))
		(negep *negep* (- negep 1)))
	       ((not (zerop (-$ (-$ one aa) one)))
		(setq *negep* (minus negep))
		(setq a aa)))
	   (setq *epsneg* a)
	   (cond ((not (or (= *ibeta* 2)
			   (zerop *irnd*)))
		  (setq a (//$ (*$ a (+$ one a))
			       (+$ one one)))
		  (cond ((not (zerop (-$ (-$ one a) one)))
			 (setq *epsneg* a)))))
	   (setq *machep* (- (minus *it*) 3)))
      (do ((aa a (*$ aa beta))
	   (machep *machep*
		   (+ machep 1)))
	  ((not (zerop (-$ (+$ one aa) one)))
	   (setq a aa)
	   (setq *machep* machep
		 *eps* aa)))
      (cond ((not (or (= *ibeta* 2)
		      (zerop *irnd*)))
	     (setq a (//$ (*$ a (+$ one a))
			  (+$ one one)))
	     (cond ((not (zerop (-$ (+$ one a) one)))
		    (setq *eps* a)))))
      (setq *ngrd* 0)
      (cond ((and (zerop *irnd*)
		  (not (zerop (-$ (*$ (+$ one *eps*) 1.0) 1.0))))
	     (setq *ngrd* 1)))
;;		  (let ((betainsq (* betain betain)))           ; betain = 1/BETA
;		    (do ((i 0)                                  
;			 (kk 1)
;			 (zsq betainsq)
;			 (yy betain  z)
;			 (z betainsq zsq)
;			 (aa (* betainsq 1.0) (* zsq 1.0)))
;			;; check for underflow
;			((or (zerop (+ aa aa))
;			     (>= (abs z)
;				 yy))
;			 (setq y yy
;			       k kk
;			       a aa)
;			 (cond ((= *ibeta* 10.)        ; for decimal machines 
;				(do ((iz *ibeta* (* iz *ibeta*))
;				     (iexp 2 (1+ iexp)))
;				    ((< k iz)
;				     (setq *iexp* iexp)
;				     (setq mx (- (+ iz iz) 1)))))
;			       (t (setq *iexp* (1+ i)
;					mx (+ k k))))) 
;		      (setq i (1+ i)
;			    kk (+ kk kk)
;			    zsq (* z z))))

      (let ((betainsq (*$ betain betain)))
	   (do ((i 0)
		(kk 1)
		(zsq betainsq)
		(z betainsq zsq))
;		(z betainsq (*$ z z))
		(aa (*$ betainsq one) (*$ zsq one))
		(yy betain z))
	       ((or (zerop (+$ aa aa))
		    (greaterp (abs z)
			      yy))
		(setq y yy)
		(setq k kk)
		(setq a aa)
		(cond ((= *ibeta* 10.)
		       (do ((iz *ibeta* (* iz *ibeta*))
			    (iexp 2 (1+ iexp)))
			   ((< k iz)
			    (setq *iexp* iexp)
			    (setq mx (- (+ iz iz) 1)))))
		      (t (setq *iexp* (+ i 1)
			       mx (+ k k))))) 
	       (setq i (1+ i)
		     kk (+ kk kk)
		     zsq (* z z)))) ; pw
      (do ((xmin y yy)
	   (yy y (*$ yy betain))
	   (aa (*$ y one) (*$ yy one))
	   (kk k (+ kk 1)))
	  ((or (zerop (+$ aa aa))
	       (not (lessp (abs yy)
			   xmin)))
	   (setq *xmin* xmin)
	   (setq k kk)
	   (setq a aa)
	   (setq *minexp* (minus k))
	   (cond ((not (or (greaterp mx (- (+ k k) 3))
			   (= *ibeta* 10.)))
		  (setq mx (+ mx mx))
		  (setq *iexp* (1+ *iexp*))))
	   (setq *maxexp* (+ mx *minexp*))
	   (setq i (+ *maxexp* *minexp*))
	   (cond ((and (= *ibeta* 2)
		       (zerop i))
		  (setq *maxexp* (1- *maxexp*))))
	   (cond ((greaterp i 20.) (setq *maxexp* (1- *maxexp*))))
	   (cond ((not (= a y))
		  (setq *maxexp* (- *maxexp* 2))))
	   (setq *xmax* (-$ one *epsneg*))
	   (cond ((not (= (*$ *xmax* one) *xmax*))
		  (setq *xmax* (-$ one (*$ beta *epsneg*)))))
	   (setq *xmax* (//$ *xmax* (*$ (*$ (*$ beta beta) beta) *xmin*))) 
	   (setq i (+ (+ *maxexp* *minexp*) 3))
	   (cond ((greaterp i 0)
		  (do ((j i (1- j)))
		      ((zerop j))
		      (cond ((= *ibeta* 2)
			     (setq *xmax* (+$ *xmax* *xmax*)))
			    (t (setq *xmax* (*$ *xmax* beta)))))))))))